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FIGURE 2.1 Stages of a bacterial growth curve. 




2.1 GROWTH MODELS 

2.1.1 Introduction 

The concept of the primary model is fundamental to the field of predictive micro- 
biology (see the definition of a model in the Preface). A primary model for microbial 
growth aims to describe the kinetics of the process with as few parameters as 
possible, while still being able to accurately define the distinct stages of growth. A 
typical bacterial growth curve is shown in Figure 2.1. When the increase in popu- 
lation density (usually defined as the base 10 logarithm of cell numbers) is plotted 
against time, the resulting curve usually has four phases, referred to respectively as 
the lag, exponential, stationary, and death or decline phases. 

In the only book published thus far that is devoted exclusively to the field of 
predictive microbiology, McMeekin et al. 1 provide an excellent review and discus- 
sion of the classical sigmoid growth functions, especially the modified logistic and 
Gompertz equations. As they point out, these are empirical applications of the 
original logistic and Gompertz functions. They lack mechanistic interpretability 
though the original logistic and Gompertz functions are considered mechanistic 
models. Over the last decade, a new generation of bacterial growth curve models 
have been developed that are purported to have a mechanistic basis: for example, 
the Baranyi model, 23 the Hills model, 45 the Buchanan model, 6 and the heterogeneous 
population model. 7 In addition to the book by McMeekin et al., other authors have 
provided reviews of microbial growth models. 38-11 

In this chapter, we will review the modified logistic and Gompertz equations as 
well as the new models that were not covered by McMeekin et al. 1 and discuss their 
applications. We will compare these models based on their performance in predictive 
microbiology applications. 
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2.1.2 The Logistic and the Gompertz Functions 

Sigmoidal functions have been the most popular ones used to fit microbial growth 
data since these functions consist of four phases, similar to the microbial growth 
curve. 9 The most commonly used are the modified logistic (Equation 2.1) and the 
modified Gompertz (Equation 2.2) introduced by Gibson et al. 12 : 



logx«) = A+ C (2.1) 



log x{t) = A + Cexp{- exp[-B(t - Af )] } (2.2) 

where x(t) is the number of cells at time t, A the asymptotic count as t decreases to 
zero, C the difference in value of the upper and lower asymptote, B the relative 
growth rate at M, and M is the time at which the absolute growth rate is maximum. 19 

The above functions use log x(t) instead of x(t) as the response variable. Thus, 
they are not simply reparameterizations of the original logistic 13 ' 14 and Gompertz 15 
functions, but are "modified" functions. The original logistic and Gompertz functions 
are considered to be mechanistic; however, the modified functions are empirical. 

The parameters of the modified Gompertz equation can be used to characterize 
bacterial growth as follows 1 : 

<? = 2.718 ■■■ 

lag time = Af-(1 //?) + - 




(2.3) 



BC/e 

exp onential growth rate = BC I e 

generation time = log(2)e / BC = 0.8183 / BC 



The expression in Equation 2.3 for lag time is different from the following Equation 
2.4 proposed by Gibson et al. 12 and other workers 16,17 : 



lag time = M (2.4) 

B 



As explained by McMeekin et al., 1 Equation 2.3 is a more general and correct 
expression for the lag time. 

In order to simplify the fitting process, reparameterized versions of the Gompertz 
equation have been proposed 18,19 : 
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where A = log 10 -x (log 10 cfu x ml [ ), x is the initial cell number, C the asymptotic 



increase in population density (log 10 cfu x ml *), R g the growth rate (log 10 cfu h '), 
and X is the lag -phase duration (h). 



2.1.2.1 Applications of the Logistic Model 

There have been limited examples of fitting of microbial growth data using the 
logistic function, since the Gompertz function, which is asymmetric about the point 
of inflection unlike the logistic function, 9 ' 20 ' 21 is generally preferred. Some recent 
examples include modeling of fish spoilage 22-24 and colony diameter of fungi. 25 A 
variation of the logistic model with a breakpoint at the transition between the lag 
phase and the exponential phase has also been used to model the lag phase of 
Listeria monocytogenes 26 




2.1.2.2 Applications of the Gompertz Equation 

The Gompertz equation has been used extensively by researchers to fit a wide variety 
of growth curves from different microorganisms. Some of the recent models devel- 
oped with the Gompertz function include those for Yersinia enterocolitica 21 Sta- 
phylococcus aureus i 2 *' 29 L. monocytogenes, 30 Vibrio parahaemolyticus 33 and Bacil- 
lus cereus. 32,33 

The Gompertz function has also been applied to growth curves based on turbidity 
data 34 ; mixed cultures of Pseudomonas spp. and Listeria spp. 35 ' 36 ; Lactobacillus 
curvatus 31 ', spoilage of vegetables, 38 beer, 39 and meat 40 ; and germination and growth 
of Clostridium botulinum.. A{ 

There are, however, some limitations associated with the use of the Gompertz 
function. The Gompertz rate (|i max ) is always the maximum rate and occurs at an 
arbitrary point of inflection 42-44 ; thus the generation time can be underestimated by 
as much as 13%. 31 In addition, since the slope of the function cannot be zero, the 
lower asymptote must be lower than the inoculum level, giving a negative X for 
some data sets. 43 Another limitation is that, in order to get a good fit, experimental 
data are required over the whole growth range. 121 




2.1.3 Baranyi Model 

In a series of papers, 2 ' 3,10 Baranyi and coworkers introduced a mechanistic model 
for bacterial growth. Briefly, the lag phase is attributed to the need to synthesize an 
unknown substrate q that is critical for growth. Once cells have adjusted to the new 
environment, they grow exponentially until limited by restrictions dictated by the 
growth medium; thus: 
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FIGURE 2.2 Example of a growth curve generated by the Baranyi and McKellar models. 
Parameters are defined in the text. 



dx 



q(t) 



r 




dt q(t) + 1 



•H 



max 



1- 



V 



' x(t) A 



\ max / 



m\ 



X(t) 



(2.6) 



/ 




where x is the number of cells at time t, x max the maximum cell density, and q(t) is 
the concentration of limiting substrate, which changes with time: 



dq 
dt 



= H 



max 



q(t) 



(2.7) 



The initial value of q (q ) is a measure of the initial physiological state of the cells. 
A more stable transformation of q may be defined as: 
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oy 



The parameter m characterizes the curvature before the stationary phase. When 
m-\ the function reduces to a logistic curve, a simplification of the model that is 
often assumed. Thus, the final model has four parameters: x , the initial cell number; 
h ; x max ; and |i- max . The output of this model (and the relationship between h , X, and 
|l max ) is shown in Figure 2.2. 45 

An explicit version of the Baranyi model has also been derived: 



2004 by Robin C. McKellar and Xuewen Lu 













~V 



1237_C02.fm Page 26 Wednesday, November 12, 2003 12:34 PM 






y(t) = y Q + li max A(t)--ln 

m 



f 



m\lA(t) 



1 + 



-1 



\ 



M>' ma x-yo) 



V 



/ 



(2.9) 




A(t) = t + — In 
v 



e + q 



!0 



V 



l + # 



(2.10) 



o 7 



where y(t) = In jc(f) , y = In x , and v is the rate of increase of the limiting 
substrate, generally assumed to be equal to \i 



max 



2.1.3.1 Applications of the Baranyi Model 

Since its inception in the early 1990s, the Baranyi model has been used extensively 
to model microbial growth. The popularity of this model has been facilitated by the 
availability of two programs: DMFit, an Excel add-in; and Micro Fit, a stand-alone 
fitting program, distributed by the Institute of Food Research in the U.K. 
(http://www.ifr.bbsrc.ac.uk/Safety/DMFit/default.html). The model was used for 
growth modeling of a wide variety of microorganisms, the results of which are 
included in the Food MicroModel software. 46 Some recent applications were related 
to Listeria monocytogenes, 47,48 B. cereus 49 Escherichia coli, 50 Y. enterocolitica, 51 
increasing colony diameter of heat-resistant fungi, 52 and spoilage in green asparagus 
and vegetable salad. 53 ' 54 

One of the advantages of the Baranyi model is that it is readily available as a 
series of differential equations that allow modeling in a dynamic environment, 
generally resulting from nonisothermal temperature profiles. This form of the model 
was used to describe the behavior of E. coli at suboptimal temperatures, 55 and to 
develop and validate a dynamic growth model for L. monocytogenes in fluid whole 
milk. 5657 It has also been used to study the influence of either slowly 58 or rapidly 59 
changing temperature on the growth of L. monocytogenes and Salmonella. 




2.1.4 Hills Model 

A general theory of spatially dependent bacterial growth in heterogeneous systems 
was developed by Hills and coworkers. 45 This was achieved by combining a struc- 
tured-cell kinetic model with reaction-diffusion equations describing transport of 
nutrients. 4 The model was based in part on the concept of DNA synthesis and cell 
division being dependent on the excess cell biomass. 

Assume M is the total biomass in the culture and N is the total number of cells 
in the culture. It can be shown that for inoculation with stationary-phase cells, 



Af(*) = Af(0)exp(Af) 

N(t) = N(0)[k n exp(Ar) + Aexp(-£/)] / (A + k n ) 



(2.11) 





A and k are rate constants; in general, they depend on all the environment factors. 
The expression for N(t) in Equation 2.1 1 has a much simpler form than the empirical 
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Gompertz function for fitting population growth, being a biexponential function 
where the second term, involving the rate of DNA synthesis, gives rise to the 
observed lag behavior. The lag time and the doubling time have the following 
relationships: 

r LAO =A- 1 log[l + (A/^)] 

(2.12) 
t LAG /t D =(ln2)- l log[l + (A/k n )] 

This shows that if the rate constants A and k n have similar activation energies, the 
ratio of lag to doubling time should be nearly independent of temperature. This 
model takes no account of possible lag behavior in the total biomass (M). 

The above model can also be generalized to spatially inhomogeneous systems 
such as food surfaces. 4 If more detailed kinetic information on cell composition is 
available, more complex multicompartment kinetic schemes can be incorporated. A 
two-compartment kinetic model of bacterial population dynamics has been devel- 
oped that is capable of describing the phenomena of lethal and sublethal injury, 
resuscitation, and transient conditions. A more general three-compartment kinetic 
model has been developed to interpret lag behavior in total biomass. These models 
can be further generalized to describe growth in spatially heterogeneous systems. 5 

2.1.4.1 Applications of Hills Model 

There have been few applications of the Hills model. The above two-compartment 
kinetic cell model was shown to fit batch-growth data for L. monocytogenes 4 and 
for Salmonella typhimurium. 5 More recently, the model was used for modeling viable 
counts of S. typhimurium in gel cassettes. 60 

2.1.5 Buchanan Three-Phase Linear Model 

Buchanan et al. 6 proposed a three-phase linear model. It can be described by three 
phases: lag phase; exponential growth phase; and stationary phase: 

Lag Phase: 

For*<* MG , N t = N 

Exponential Growth Phase: 

For t LAG < t < t MAX , N t = N + \L(t - t LAG ) (2. 1 3) 

Stationary Phase: 

For t > t MAX , N t = N MAX 

where N t is the log of the population density at time t (log cfu ml -1 ); JV the log of 
the initial population density (log cfu ml -1 ); N MAX the log of the maximum population 
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FIGURE 2.3 Example of a growth curve generated by the Buchanan model. Parameters are 
defined in the text. 




density supported by the environment (log cfu ml -1 ); t the elapsed time; t LAG the 
time when the lag phase ends (h); t MAX the time when the maximum population 
density is reached (h); and |i is the specific growth rate (log cfu ml -1 Ir 1 ). The three- 
phase model is illustrated in Figure 2.3. 

In this model, the growth rate was always at maximum between the end of the 
lag phase and the start of the stationary phase. The |i was set to zero during both 
the lag and stationary phases. The lag was divided into two periods: a period for 
adaptation to the new environment (t a ) and the time for generation of energy to 
produce biological components needed for cell replication (t m ). Thus, the lag phase 
is given by: 




LAG a m 



(2.14) 



This implies that t a and t m can be estimated from data fitted with the linear model 
using the following relationships 6 : 



t = generation time 

t a = t LAG - generation time 



(2.15) 



2.1.5.1 Applications of the Buchanan Model 

Surprisingly, this simple model has not been used extensively for fitting growth data. 
The original authors used the three-phase version of the model to fit experimental 
data for E. coli 0157:H7. 6 As there is often little interest in modeling the stationary 
phase, a modified two-phase version has been proposed that fits only the lag and 
exponential phases. In a series of papers published in 1999, Oscar used the two- 
phase model to fit growth data for S. typhimurium in brain heart infusion broth, 61 
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and on cooked chicken 62 and ground chicken 63 breast meat. Fitting was accomplished 
using a useful nonlinear regression software package called Prism (GraphPad Soft- 
ware, San Diego, CA) in which an if-then statement defines the model: 

N t =N o +mt^t IAa ,0,\i-(t-t IAa )] (2.16) 

with symbols defined in Equation 2.13. A two-phase model was also used to model 
growth of E. coli 0157:H7. 64 

2.1.6 McKellar Model 

One of the limitations of existing models is that they all assume a homogeneous 
population of cells. A heterogeneous population model was recently proposed in 
which growth was expressed as a function of two distinct cell populations. 7 Cells 
can exist in one of two "compartments" or states: growing or nongrowing. All growth 
was assumed to originate from a small fraction of the total population of cells that 
are present in the growing compartment at t = 0. Subsequent growth is based on the 
following logistic equation: 



dG ., 

= G-\L 

dt 






N 
V v MAxy 



(2.17) 




where G is the number of growing cells in the growing compartment. The majority 
of cells were considered not to contribute to growth, and remained in the nongrowing 
compartment, but were included in the total population. While this is an empirical 
model, it does account for the observation that growth in liquid culture is dominated 
by the first cells to begin growth, and that any cells that subsequently adapt to growth 
are of minimal importance. 7 

This model has an interesting relationship with the Baranyi model. It is derived 
from a different initial premise, that microbial populations are heterogeneous rather 
than homogeneous. It is based on two populations of cells that behave differently, 
rather than a single population. The sum of the two populations effectively describes 
the transition from lag to exponential phase, and defines a new parameter G , the 
initial population capable of growth. Reparameterization of the model led to the 
finding that a relationship existed between |i. max and X, which is shown in Figure 
2.2, 7 and which had been derived by Baranyi from a more mathematical argument. 3 
Baranyi 65 later supported the geometric relationship in Figure 2.2, and stated that 
the initial physiological state of the whole population could reside in a small sub- 
population. Thus, the McKellar model constitutes a simplified version of the Baranyi 
model, and has the same parameters. 

The concept of heterogeneity in cell populations was extended further to the 
development of a combined discrete-continuous simulation model for microbial 
growth. 66 At the start of a growth simulation, all of the cells were assigned to the 
nongrowing compartment. A distribution of individual cell lag times was used to 
generate a series of discrete events in which each cell was transferred from the 
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nongrowing to the growing compartment at a time corresponding to the lag time for 
that cell. Once in the growing compartment, cells start growing immediately accord- 
ing to Equation 2.17. The combination of the discrete step with the continuous 
growth function accurately described the transition from lag to exponential phase. 
This model was further modified to include a continuous adaptation phase prior to 
the discrete event. 67 A new physiological state parameter p was proposed that 
represents the mean of the initial individual cell physiological states. This model is 
dynamic in both the lag and exponential phases, and thus is useful for simulating 
the behavior of individual cells in a changing environment. 

2.1.6.1 Applications of the McKellar Model 

This model has not been used extensively for modeling microbial growth partly 
because of its similarity to the Baranyi model. It is also a compartmental model, 
and as such cannot be fitted easily using conventional nonlinear regression programs. 
This model was fitted to data for growth of L. monocytogenes at 5 to 35°C, and 
compared to the Gompertz model. 7 Values for |l max were slightly higher with this 
model, and X were generally shorter than found with the Gompertz model. Goodness- 
of-flt analysis suggested that the McKellar model generally fit the data better than 
the Gompertz. 

2.1.7 Other Models 

There have been a large number of alternative models proposed for modeling micro- 
bial growth. Many of the earlier ones have been thoroughly discussed by McMeekin 
et al., 1 and will not be discussed further here. 

Whiting and Cygnarowicz-Provost 68 suggested a quantitative four-parameter 
model for the germination, growth, and decline of C. botulinum, and the growth of 
L. monocytogenes. Jones and Walker 69 developed an equation to predict growth, 
survival, and death of microorganisms based on data obtained using Y. enterocolitica 
in varying pH and sodium chloride concentrations at different temperatures. Van 
Impe et al. 70 proposed a dynamic first-order differential equation to predict both 
microbial growth and inactivation, with respect to both time and temperature. We 
are expecting more accurate and more mechanistic primary models when people 
gain more knowledge on the kinetics of individual cells and behavior of bacteria. 
Recently, a series of three models has been proposed in which |l can increase, remain 
constant, or decrease with time. 71 The latter two models bear some resemblance to 
those discussed earlier; however, the concept of (J increasing with time was designed 
to accommodate the observation that recombinant E. coli initially grew rapidly in a 
bioreactor because of high substrate concentrations. 

2.1.8 Examples of Growth Model Fitting 

It seems appropriate at this point to provide an example of how some of the more 
popular and useful functions may be used to fit experimental growth data. The data 
selected (taken from an earlier study 7 ) were for the growth of L. monocytogenes at 
5°C (Table 2.1). The models used in this comparison were Gompertz using Equation 
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TABLE 2.1 

Growth Data for Listeria monocytogenes 

at 5°C 



Time (d) 



6 

24 

30 

48 

54 

72 

78 

99 
126 
144 
150 
168 
174 
191 
198 
216 
239 
266 
291 
316 
336 
342 
360 
384 



log cfu ml 1 



4.8 
4.7 
4.7 
4.7 
4.9 
5.1 
5.3 
5.4 
5.9 
6.3 
6.9 
6.9 
7.2 
7.3 
7.7 
7.8 
8.3 
8.8 
9.1 
9.2 
9.3 
9.7 
9.7 
9.7 
9.5 






2.5, Baranyi using Equation 2.6 and Equation 2.7, McKellar using Equation 2.17, 
and Buchanan using Equation 2.13. Nonlinear regression analysis was done using 
the ModelMaker® software (Modelkinetix, Old Beaconsfield, Bucks, U.K., 
www.modelkinetix.com), which uses the Runge-Kutta method for solving differen- 
tial equations. Initial parameter estimates were made using the simplex method, and 
regression was performed using the Marquardt algorithm. The Baranyi and McKellar 
models gave values for |i max directly, since they were in the form of differential 
equations, and modeled the cell number rather than log 10 cfu ml -1 . The Gompertz 
and Buchanan models were applied directly to log 10 cfu ml -1 data, and thus the rate 
parameter (/? ) obtained from the fitting had to be converted to |i max using the 
relationship \l max = R -In 10. The X parameter for the Gompertz and Buchanan 
models was obtained directly from the fitting, while the values for the Baranyi and 
McKellar models were derived from the h parameter values using the following 
relationship: h = \i max • X . The Baranyi model (Baranyi MF ) was also fitted using the 
MicroFit software, in which the model was reparameterized to fit X directly. The 
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TABLE 2.2 

Results of Model Fitting to Growth Data 

Model 3 \i (h" 1 ) k (d) Log x Q Log x max DF RMSE 



Baranyi MF 


0.050 


46.9 


4.68 


9.57 


21 


0.100 


McKellar 


0.049 


44.9 


4.63 


9.57 


21 


0.112 


Gompertz pz 


0.054 


54.7 


4.68 


10.0 


21 


0.119 


Gompertz 


0.058 


68.4 


4.76 


9.87 


21 


0.139 


Buchanan 


0.048 


53.8 


4.84 


9.49 


20 


0.157 


Bar any i 


0.056 


61.4 


4.72 


9.32 


21 


0.179 



Note: DF = residual degrees of freedom; RMSE = root mean square error. 

a The McKellar, Gompertz, Buchanan, and Baranyi models were fit using the 
ModelMaker software. Baranyi MF is the Baranyi model fit using the MicroFit 
software. Gompertz pz is the Gompertz model fit using the Prism software. 



Gompertz model (Gompertz PZ ) was also fitted using Prism™ Version 3.03 (GraphPad 
Software, San Diego, CA, www.graphpad.com). 

The results of the various fitting approaches are given in Table 2.2. The root 
mean square error (RMSE) was taken as the measure of goodness of fit, as suggested 
by Ratkowsky (Chapter 4). The models are placed in order of increasing RMSE. 

The best model was Baranyi MF , with the lowest RMSE. In contrast, the poorest 
fit was with the Baranyi model using ModelMaker, which also gave a higher |l max 
and X than did the Baranyi MF model. The next best model was the McKellar, with 
parameter values close to those for Baranyi MP The Gompertz model fitted using 
either Prism or ModelMaker gave larger |l max and X values than did the Baranyi MF 
and McKellar models, and the highest values of log x max among all other models. 
The Buchanan model gave the lowest value of |l max of all the models, and a shorter 
A, than all except the Baranyi MF and McKellar models. 

The output of the four models fitted with ModelMaker is also shown in Figure 
2.4. The steeper slope (|i max ) of the Gompertz and Baranyi models may be observed. 
The greatest difference between models occurred during the late-log early-stationary 
phase. The Gompertz model never reaches a plateau, which reflects its higher log, 
x max (Table 2.2). As expected, the Buchanan model has a sharp breakpoint, while 
the transition to the stationary phase appears smoother with both the McKellar and 
Baranyi models. 

The results of the nonlinear regression fitting described above emphasize an 
important point: there is no single solution for nonlinear regression, in contrast to 
linear regression. The iterative approach used in nonlinear regression is dependent 
on the parameter starting values, and may find local, rather than global, optimum 
values. In addition, different software packages use different procedures for fitting, 
and thus the results obtained (such as those above) should be considered comparative 
rather than absolute. The fitting results do show that, while there are differences 
between the models and the software used, the parameter differences are often slight. 
It is worth noting that estimates of X range from 44.9 to 68.4 days. Further discussion 
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FIGURE 2.4 Comparison of growth models fitted to viable count data of Listeria monocy- 
togenes grown at 5°C. 

on the difficulties in modeling X, and the role of the physiological state, can be found 
in Chapter 9, and a more complete discussion of model fitting can be found in 
Chapter 4. 

2.1.9 Comparison of Existing Models 

Zwietering et al. 18 statistically compared several different modified sigmoidal func- 
tions (Logistic, Gompertz, Richards, Schnute, and Stannard) using the Mest and the 
F-test. In most of the cases, the modified Gompertz expression was regarded as the 
best model to describe the growth data both in terms of statistical accuracy and ease 
of use when compared to other sigmoidal functions. 

Baranyi et al. 2 compared the output of their model with that of the Gompertz, 
and concluded that the goodness of fit was generally at least as good. They also 
showed that their model gave estimates for lag and growth rate that were slightly 
lower than in the Gompertz case. Baranyi et al. also compared their model to those 
of Hills 10 and Buchanan 72 and stated that these models are special cases of the 
Baranyi model. Baranyi argues that the Buchanan model has merit in its simplicity, 
but that the model lacks the capability of simulating dynamic behavior. 72 Buchanan 
et al. 6 asserted that their three-phase model is comparable to, and more robust than, 
either the Gompertz or the Baranyi models, especially when experimental data were 
minimal. The three-phase linear and Baranyi models predicted similar maximum 
population densities. These values were typically smaller than the values provided 
by the Gompertz model. Garthright 44 strongly supports the three-phase model, and 
points out its superiority in describing the lag and exponential phases as compared 
to the Gompertz. He concludes that the nonlinear approach does not achieve any 
advantage over the three-phase linear approach for environmental applications. This 
model appears particularly appropriate for modeling conditions where growth is 
poor, and an upper asymptote cannot be accurately fixed. The Baranyi model and 
the McKellar model can also be used when stationary-phase data are lacking. 

2004 by Robin C. McKellar and Xuewen Lu 











1237_C02.fm Page 34 Wednesday, November 12, 2003 12:34 PM 







Other comparisons between growth models have been made. A comparison of 
the logistic, Gompertz, and Baranyi models for fish spoilage showed that the logistic 
function was similar to the Baranyi but easier to fit. 73 A comparison between Gomp- 
ertz and Baranyi models gave better fit with the Baranyi model, and a higher growth 
rate with Gompertz. 74 The Gompertz function was found to be more appropriate 
than the Baranyi model for monitoring C0 2 evolution as an indicator of bacterial 
growth. 75 Other workers have compared the Baranyi and Gompertz models, and have 
concluded that the Baranyi function gave better parameter estimates as compared to 
the Gompertz. 76 

At the present time it is not possible to select one growth model as the most 
appropriate representation of bacterial growth. If simple is better, then the three- 
phase model is probably sufficient to represent fundamental growth parameters 
accurately. 4477 There does appear to be general agreement in relationship to under- 
lying principles, and emphasis should be placed on the development and use of 
models and parameters that can be easily understood by food microbiologists. 77 
However, in spite of Garthright's assertion that straight line simplicity is sufficient 
to model growth, 44 the development of more complex models (and subsequently 
more mechanistic models) will depend on an improved understanding of cell behav- 
ior at the physiological level. 




2.2 SURVIVAL MODELS 

2.2.1 Introduction 

Our ability to understand and model the survival of pathogens in foods or during 
processing of food is critical to the safety of the food supply. Thus, models to 
describe microbial death due to heating have been used since the 1920s, and con- 
stitute one of the earliest forms of predictive microbiology. Much of the early work 
centered around the need to achieve destruction of C. botulinum spores in low acid 
canned foods, 17 ' 78,79 and much effort has been put towards characterizing the kinetics 
of spore inactivation. In this section of the chapter we will focus on the evolution 
of survival modeling from the classical linear approach to the more complex models 
required to describe inactivation curves that deviate from linearity. 

2.2.2 Classical Linear Models 

It has always been assumed that spore inactivation follows simple first-order reaction 
kinetics under isothermal conditions: 



dS t 

—^ = -k'S t (2.18) 

dt 

where S t is the survival ratio (NJN Q ) and k' is the rate constant. Thus the number of 
surviving cells decreases exponentially: 
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S,=e 



-k't 



(2.19) 



and when expressed as log 10 , gives: 



log S. = -kt 



(2.20) 



where k = k'/ln 10. The well-known D- value (time required for a 1-log reduction) 
is thus equal to Ilk, where k is the slope (Figure 2.5). The D-values can also be 
expressed as: 



D -value = 



log N, - log N t 



(2.21) 



When log D-values are plotted against the corresponding temperatures, the reciprocal 
of the slope is equal to the z-value, which is the increase in temperature required 
for a 1-log decrease in D-value (Figure 2.5; inset). The rate constant can also be 
related to the temperature by the Arrhenius equation: 



k = N e 



F A 
RT 



(2.22) 




where E a is the activation energy, R the universal gas constant, and T is the temper- 
ature in Kelvin. 

From the first-order reaction it is not possible to achieve complete destruction 
of all C. botulinum spores in a given volume of product; one spore will always be 





Time (min) 

FIGURE 2.5 Geometric description of D- and z-values. (From McKellar, R.C., Modelling 
the effectiveness of pasteurization, in Dairy Processing: Maximizing Quality, Smit, G., Ed., 
CRC Press Inc./Woodhead Publishing, 2003 pp. 104-129. With permission.) 
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left in a can if a sufficient number of cans are examined. Thus it is generally assumed 
that a 12-log reduction (also known as 12D) is sufficient to achieve "commercial 
sterility," or an acceptable level of risk of survival of C. botulinum. Knowledge of 
the D-values of representative strains allows the determination of the F -value, which 
is the time required to achieve 12D, assuming a z- value of 10°C. At 121°C, F is 
equal to 2.5 min for most strains of C. botulinum. 11 

Comparable standards for other food-borne pathogens do not exist; however, it 
is generally accepted that a 4- or 5 -log reduction is considered adequate, depending 
on the product. An extensive amount of work has gone into the determination of D- 
and z-values for various pathogens. Thermal stability of pathogens such as L. mono- 
cytogenes* salmonellae, 81 and E. coli 0157:H7 82 has been well documented and 
summarized in recent reviews. 

2.2.3 Nonlinear Models 
2.2.3.1 Nonlinearity Issues 

The canning industry has enjoyed an enviable record of safety, and thus the concept 
of logarithmic death of microorganisms has persisted, and is now considered 
accepted dogma. In spite of this, nonlinear survival curves were reported for some 
bacteria almost 100 years ago. 83 In general there are two classes of nonlinear curves; 
those with a "shoulder" or lag prior to inactivation, and those that exhibit tailing. 
These two phenomena may be present together, or with other observed kinetics such 
as biphasic inactivation. A wide variety of complex inactivation kinetics have been 
reported, and several of these are shown in Figure 2.6. The theoretical basis for 
assuming logarithmic behavior for bacteria is based on the assumptions that bacterial 
populations are homogeneous with respect to thermal tolerance, and that inactivation 
is due to a single critical site per cell. 83 Both of these assumptions have been 
questioned, and thus concerns have been raised regarding the validity of extrapolation 
of linear inactivation curves. 84 - 85 

Stringer et al. 82 have summarized the possible explanations for nonlinear kinetics 
into two classes: those due to artifacts and limitations in experimental procedure 
and those due to normal features of the inactivation process. The first class encom- 
passes such limitations as variability in heating procedure; use of mixed cultures or 
populations; clumping; protective effect of dead cells; method of enumeration; and 
poor statistical design. The second class includes such situations as possible multiple 
hit mechanisms; natural distribution of heat sensitivity; and heat adaptation. These 
two classes roughly parallel the two concepts reviewed by Cerf 85 to explain tailing 
in bacterial survival curves. The first of these (the "mechanistic" approach) also 
makes the assumption of homogeneity of cell resistance and proposes that thermal 
destruction follows a process analogous to a chemical reaction. In this approach, 
deviations from linearity are attributed mainly to artifacts; however, tailing is also 
related to the mechanism of inactivation or resistance. In the "vitalistic" approach, 
it is assumed that the cells possess a normal heterogeneity of heat resistance; thus 
survival curves should be sigmoidal or concave upward. 85 
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FIGURE 2.6 Examples of thermal death curves: (a) lag or shoulder, with either linear (dotted 
line), power law where p > 1 (broken line), or monophasic logistic (solid line) models; (b) 
concave with power law where/? < 1; (c) biphasic logistic; and (d) sigmoidal. 

There has been considerable controversy between the two schools of thought, and 
the literature is divided on the validity of nonlinear survival curves as representing 
the true state of the cell population. There is certainly evidence that inconsistencies 
in experimental protocols or the use of incorrect media can lead to artifacts; however, 
there is little convincing evidence that clumping of cells or the protective effect of 
dead cells is consistently responsible for nonlinear survivor curves. The current belief 
is, notwithstanding some contribution by artifacts, that cells do exhibit heterogeneity 
in thermal sensitivity, and the majority of modeling approaches now make this assump- 
tion. There is also inconsistency in actually defining what is meant by an artifact. If 
one assumes that an artifact in this context is anything that interferes with obtaining 
a linear death curve, then many of the situations currently classified as artifacts may 
be natural behavior of cell populations. This is particularly obvious in the study of 
spore inactivation where standardized suspensions are difficult to obtain, and much 
effort has been expended to remove artifacts such as genetic variants. The difficulty 
in obtaining linear kinetics may be a signal that, in most cases, nonlinearity is the norm. 

The current theories of microbial inactivation must be revisited in light of recent 
improved understanding of the effect of heat on microorganisms. We now know that 
cells do not exist simply as alive or dead, but may also experience various degrees 
of injury or sublethal damage, which may give rise to apparent nonlinear survival 
curves. 82 The induction of heat resistance in food-borne pathogens due to expression 
of heat shock proteins has been extensively documented in recent years, and may 
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contribute to apparent nonlinearity, particularly tailing. 8286 ' 87 Thus it appears impor- 
tant to model the actual conditions or situations experienced by bacteria in foods 
rather than relying on simplifications. Survival modeling should also include a more 
complete understanding of the molecular events underpinning microbial resistance 
to the environment. 

It seems likely that heterogeneity within bacterial populations is responsible in 
most cases for nonlinear survival curves, and most recent attempts to model survival 
employ distributions. The use of distributions to account for nonlinearity is not new; 
log normal distributions had been suggested for this purpose as early as 1942. 83 
Other distributions such as logistic, gamma, and Weibull have also been suggested; 
Weibull is the favored approach at the moment (see later). There is no complete 
agreement on the use of distributions, 83 and it is clear that this approach cannot 
adequately account for changes in heat resistance occurring during heating. 

Our lack of understanding of the key physiological aspects of microbial inactivation 
and the complexities of nonlinear behavior suggest that a truly mechanistic model for 
thermal inactivation will not be developed in the near future. One approach to quanti- 
tating bacterial survival might be the thermal death point concept common to the 
canning industry. This approach allows one to define the conditions required to achieve 
a target log reduction, and makes no statement regarding the kinetics of that destruction. 
This approach has a number of attractive advantages; however, it would still be influ- 
enced by such artifacts as changes in heat resistance of a culture and cell injury. 83 

2.2.3.2 Shoulder/Tail Models 

2.23.2.1 Linear Approach 

Inactivation curves that deviate from simple exponential often have a lag or shoulder 
region prior to the exponential inactivation. This shape of inactivation curve is 
probably the most commonly experienced by researchers. A simple linear model to 
account for this behavior was developed by Whiting 88 : 




\ogN = < 



log N when < t < t L 

( n 

logAf n - — (t-t r ) when t>t f 



(2.23) 



\Dj 



where t L is the lag prior to inactivation. 

An example of the output of this model is the dotted line in Figure 2.6a. The 
advantage of this model is that linear regression can be used. This simple model has 
been used effectively to describe the nonthermal inactivation of L. monocytogenes 
as a function of organic acid and nitrite concentrations 89-92 and under reduced 
oxygen. 93 A similar two-phase linear model was described for thermal inactivation 
of L. monocytogenes by Breand et al. 94 

It is quite common for the lag or shoulder region of the survival curve to be 
highly variable. This makes it difficult to develop secondary models to describe the 
influence of the environment on the lag. Thus, survival using this model is often 
described as the time required for a 4-log reduction (r 4D ) 89 ' 92 ' 95 : 
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t 4D= t L +4 ' D 



(2.24) 



2.2.3.2.2 Nonlinear Approach 

Complex inactivation kinetics requires the use of nonlinear functions. It should be 
noted that nonlinearity as it relates to mathematical functions means that the param- 
eters in the equation are nonlinear; the resulting curve may or may not appear linear. 
Linear regression can be easily performed by most spreadsheet programs; however, 
nonlinear regression is an iterative process that is supported by more specialized 
software. These software packages are readily available; thus considerable advances 
have been made in the development of nonlinear models. 

Another of the more common shapes of survival curves is the concave curve, 
which has no lag, and a single, tailing population (Figure 2.6b). This function is 
best represented by the power law: 



log 






D 



(2.25) 




where p is the power. A concave curve is produced when p < 1 (Figure 2.6b), and 
a convex (or shoulder) shape results from p > 1 (broken line in Figure 2.6a). A 
power law function has been used to model curvature in survival curves for Entero- 
coccus faecium 96 and alkaline phosphatase 97 in milk. Other, seemingly novel, func- 
tions that have been derived to fit concave survival curves are really in fact power 
law functions. 98,99 

Tailing survival curves can also be represented by the exponentially damped 
polynomial model. In this model, deviation from simple linear kinetics, experienced 
while heating Staphylococcus aureus in skim milk, was fitted with the nonlinear 
function 100 : 




N 
log — = -kte~ h 



(2.26) 



where k is the rate coefficient and X is the damping coefficient. 

As discussed earlier, a logistic equation may be used in growth modeling to 
modify the simple exponential growth to account for limiting the maximum popu- 
lation size as a result of nutrient limitation. In the same way, a logistic function can 
be used to account for death being limited by the amount of some stress factor or 
damage to the cell. 101 This "mirror image" of the logistic function is called the Fermi 
equation, and is used for sigmoidal decay curves, which are symmetric: 



i N i 



TV, 



o 



1 i ~ bt i 



1 + e 



b(t-t L ) 



(2.27) 



where N is the population (cfu ml l ) surviving at time t\ N is the population surviving 
at time 0; b is the maximum specific death rate; and t L is the lag phase prior to 





2004 by Robin C. McKellar and Xuewen Lu 











~V 



1237_C02.fm Page 40 Wednesday, November 12, 2003 12:34 PM 





inactivation. This equation has been modified to account for situations where one 
may find both a primary, heat-sensitive population and a secondary more heat- 
resistant population 88 : 



i N i 



F(l + e" ¥i ) (1-F )(1 + <T Vi ) 



b x (t-t L ) 



d + e^'-'L') (\ + e u ^-' L >) 



b 2 (t-t L ) 



(2.28) 




where Z?j is the maximum specific death rate for the primary population and b 2 is 
the maximum specific death rate for the secondary population. Traditional D-values 
may be calculated as 2.3/b for each population. Lag phases are not always present, 
and this can be accounted for by setting the value of t L to zero. An example of the 
output of this function is given in Figure 2.6c. The biphasic logistic model has been 
used to model inactivation of spores of C. botulinum, 102 and the nonthermal inacti- 
vation of L. monocytogenes 90 - 92 ' 93 and S. aureus. 103 This model has also been applied 
to the thermal inactivation of bovine milk catalase 104 and E. faecium 96 during high- 
temperature short-time (HTST) pasteurization, and inactivation of E. faecium during 
bologna sausage cooking. 105 In situations where a single population exists, F can be 
set equal to 1 (solid line in Figure 2.6a). 

Other variations of the logistic function have been suggested. A four-parameter 
logistic model was proposed by Cole et al. 106 : 



y = a + 



co-a 




4g(t-x) 



(2.29) 



\ + e 



oj-a 



where y = log 10 cfu ml -1 ; x = log 10 time; a = upper asymptote; CO = lower asymptote; 
X = position of maximum slope; and o = maximum slope. This model was applied 
to the survival of Y. enterocolitica at suboptimal pH and temperature, 107 and the 
thermal inactivation of Salmonella typhimurium, im C. botulinum, 109 Salmonella 
enteritidis, and E. coli. no 

As was shown earlier, the asymmetric Gompertz function has considerable 
advantages when fitting bacterial growth curves. In keeping with the trend to use 
mirror images of growth functions to describe inactivation, a reparameterized form 
of the Gompertz function was suggested by Linton et al. 111 : 



log 



TV. 



= Cexp(- exp(A + Bt)) - Cexp(- exp(A)) 



(2.30) 



This function has been used to fit nonlinear survival curves of L. monocytogenes 
in buffer 111 and infant formula. 112 An example of the Gompertz function is given 
in Figure 2.6d. Other applications for the Gompertz equation include the effect of 
combined high pressure and mild heat on the inactivation of Escherichia coli and 
S. aureus in milk and poultry, 113 and the inhibition of Enterobacteriaceae and 
Clostridia during sausage curing. 114 In a similar fashion, the mirror image of the 
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TABLE 2.3 

Survival Data from Pediococcus sp. 

NRRL B2354 at 62°C 



Time (min) 


Log 


cfu ml 1 







8.4 


5 




8.3 


10 




8.1 


15 




8.1 


20 




7.7 


25 




7.3 


30 




6.9 


35 




6.4 


40 




6.2 





TABLE 2.4 

Parameter Estimates from Fitting the 
Logistic and Linear Survival Models to 
the Data in Table 2.3 




Model 


N 


Lag 


D- value 


DF 


RMSE 


Logistic 


8.11 


12.3 


12.0 


6 


0.089 


Lineai - 


8.27 


12.8 


12.7 


6 


0.102 



Baranyi growth model (see earlier) has been used for fitting nonlinear survival 
curves for the thermal inactivation of Brochothrix thermosphacta 115 and Salmonella 
enteritidis. 116 

2.2.3.2.3 Examples of Model Fitting 

It has often proven difficult to accurately fit survival data where a lag exists prior 
to inactivation. The models we have found most useful in this situation are the single- 
phase logistic (Equation 2.27) and the two-phase linear (Equation 2.23). 

These models were fitted to unpublished data on survival of Pediococcus sp. 
NRRL B2354 heated at 62°C (Table 2.3), using Prism as described above. The 
results of the fitting are shown in Table 2.4, and in Figure 2.7. The logistic model 
was slightly better than the linear, with a smaller RMSE. Because of the sharp 
breakpoint between the shoulder and exponential decay, the D-value for the linear 
model was slightly larger, while the lag phase was only marginally greater than that 
in the logistic model. As was found with growth models, there is often little to choose 
between different models; thus personal preference and experience often dictate 
which model is generally used. A more complete discussion of model fitting may 
be found in Chapter 4. 
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FIGURE 2.7 Example of fitting nonlinear survival data for Pediococcus NRRL B2354 using 
monophasic logistic (solid line) and two-phase linear (broken line) models. (From McKellar, 
R.C., Modelling the effectiveness of pasteurization, in Dairy Processing: Maximizing Quality, 
Smit, G., Ed., CRC Press Inc./Woodhead Publishing, 2003 pp. 104-129. With permission.) 




2.2.4 Distributions 

One recent development in the modeling of bacterial survival is the use of distribu- 
tions. This is based on the assumption that lethal events are probabilistic rather than 
deterministic. With a large initial population of cells, a continuous function can be 
used, much like with a chemical reaction (although a chemical reaction appears 
deterministic only because of the large number of molecules involved). The survival 
curve for a single cell is a step function, where a cell exists as either alive or dead 117 : 




S i (t) = 



1 (alive) for t < t c 
(dead) for t > t 



(2.31) 



where t c is the inactivation time. Since all cells would not be expected to die at the 
same time, values of t c would follow some sort of distribution. The Weibull distri- 
bution is used in engineering to model time to failure, and so it seems appropriate 
for modeling bacterial inactivation. The distribution of t c would then follow the 
probability density function (PDF) for the Weibull (solid line in Figure 2.8): 



PDF = 



P 
a 



r t\ 



va; 



P-i 



a 



(2.32) 



where a and P are parameters relating to the scale and shape of the distribution, 
respectively. 118 The survival curve is then the cumulative distribution function (CDF) 
(dotted line in Figure 2.8): 
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FIGURE 2.8 Probability density (solid line) and cumulative probability distribution (broken 
line) for the Weibull distribution. 



CDF = e 



a 



(2.33) 




It can be easily seen that the CDF of the Weibull distribution is essentially a 
reparameterization of the power law function (Equation 2.25). In the same fashion, 
the Fermi equation described earlier is the CDF of a log normal PDF. 21 ' 86 

The Weibull parameter (P) has a very distinct influence on the shape of the 
survivor curve. When p < 1, a concave survival curve is obtained, and when P > 1, 
the curve is convex. Interestingly, the simple exponential model described earlier is 
a special case of the Weibull distribution when p = 1, providing further support for 
the use of the Weibull distribution as an effective modeling approach for microbial 
survival. Further, the value of p can have some implications for possible mechanisms 
of inactivation. When p < 1, there is an indication that the remaining cells are more 
resistant to the treatment, while when P > 1 , an accumulation of the lethal effect is 
observed resulting in increasing rate of destruction with time. The classical D-value 
from linear survival curves can be related to the 90% percentile of the CDF, which 
is the time (t d ) required to reduce the number of microorganisms by a factor of 10 118 : 




t d =a(2.303) p 



(2.34) 



There have been a number of recent applications of the Weibull distribution to 
model survival curves for species of Bacillus and Clostridium spp., 98 Salmo- 
nella, 119 * 120 and E. coli. nx Van Boekel 118 has fitted the Weibull distribution to a large 
number of survival curves obtained from the literature. In almost all cases, the P 
values were different from 1, indicating that the classical linear model may not be 
generally applicable. Temperature had a significant effect on the a but not the p 
parameter. In order to determine if the Weibull distribution is appropriate for a 
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FIGURE 2.9 Model for spore activation and survival. 

particular survival curve, a so-called hazard plot 118121 of ln(-ln S) vs. In t should 
give a straight line. It should also be noted that, when survival curves are modeled 
using distributions, the presence of a "shoulder" can be attributed to the spread of 
a distribution being small relative to its mean or mode. 122 

2.2.5 Spores 




Modeling the inactivation of bacterial spores presents a unique problem. Spore- 
forming bacteria such as Bacillus and Clostridium spp. can exist in a dormant (spore) 
stage that is highly heat resistant. Germination of spores can be achieved by treatment 
with sublethal heat. 123 Because of the extreme heat resistance of some of these micro- 
organisms, activated spore preparations have traditionally been used to establish ster- 
ilization protocols in the canning and ultrahigh temperature industries. 124 As described 
earlier, the classical view of microbial thermal inactivation ascribes a first-order 
reaction to the process; however, it has been difficult to consistently achieve simple 
exponential inactivation with spore preparations. These variations manifest themselves 
as a shoulder on the decay curve, which has been attributed to activation of spores, 
and subsequent differences in the heat resistance of dormant and activated spores. 125 
Consistent populations of activated spores are difficult to obtain; thus the shoulder is 
often ignored, and D- values are calculated from the linear portion of the decay curve. 
More sophisticated models have been developed to account for the nonlinear 
aspects of survival curves. These include terms describing the germination of spores 
prior to inactivation (for descriptions of earlier models, see). 124-126 Figure 2.9 indi- 
cates the process of activation of dormant spores {N x ) into activated (A 2 ) spores with 
rate constant of k a . The activated spores are then inactivated by heat treatment (Ay 
at a rate equal to k d2 . The model also allows for possible inactivation of dormant 
spores (AT 3 ) at a rate equal to k d[ . All reactions are considered to be independent first- 
order. The simplest form of this model was described by Shull et al., 127 and assumes 
that only activated spores can be killed (k dl = 0) and thus: 




dN 1 
dt 



= -k A, 

a 1 



(2.35) 



dN 2 
dt 



= k A, -k n N 



(2.36) 



The model proposed by Rodriguez et al. 128 ' 129 advances the Shull model by 
assuming that the dormant spores can also be inactivated: 
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dN x 
dt 



= -(k d2 + k a )N ] 



(2.37) 



and k d2 = k dV Sapru et al. 124 further extended the model to include a different rate 
of inactivation for dormant spores (k dl ^ k d2 ). The Sapru model is more general and 
includes the Shull and Rodriguez models as special cases. This model was proposed 
for use with Bacillus stereothermophilus at sterilization temperatures, and an explicit 
form has been presented 125 : 



NAt) = NA0)e~ (ka+kdl)t 



(2.38) 



N 2 (f) = N 2 (Q)e~ kd2 ' + B-N l (0)(l- e~ m )e 



-At\ _-k d2 t 



(2.39) 



with 



B = 



A 



(2.40) 



a cl\ dl 




(2.41) 



and where N x (0) and A^ 2 (0) are the number of dormant and activated cells, respec- 
tively, at t = 0. An example output from the Sapru model is shown in Figure 2.10. 
With ^(0) at 1 x 10 8 and -/V 2 (0) at 1 x 10 5 , the initial rapid increase in surviving 
cells is the result of spore activation. This is followed by an exponential decrease 
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FIGURE 2.10 Output of model for spore activation and survival. (From McKellar, R.C., 
Modelling the effectiveness of pasteurization, in Dairy Processing: Maximizing Quality, Smit, 
G., Ed., CRC Press Inc./Woodhead Publishing, 2003 pp. 104-129. With permission.) 
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in surviving cells. This model has been expanded further to include subpopulations 
of spores having different heat resistances. 130 

2.2.6 Processing Models 

2.2.6.1 Thermal 

Thermal inactivation of microorganisms in static or batch systems is usually 
described by the D- and z-value concepts as discussed above, with temperature 
generally held constant. The situation in canning operations or continuous flow 
systems such as HTST pasteurization, sterilization, and ultrahigh temperature pro- 
cesses is somewhat more complex, due to nonisothermal conditions. In addition, the 
kinetics of inactivation in continuous systems differs from that in batch systems 
since in these systems there are additional factors such as pressure and shear forces 
that can influence microbial survival. 131 As most modern processes are continuous, 
it is necessary to have information on survival of microorganisms; however, few 
studies have been published in which laboratory or pilot plant continuous flow 
systems have been studied. 131 

In order to deal with nonisothermal conditions, Bigelow's 132 model has been the 
standard for the low-acid canned food industry for many decades. In this approach, 
the processing time F is determined by integrating the exposure time at various 
temperatures, T[f], to time at a reference temperature, T Ref 133 : 





J 



(T(Q-T Ref ) 

F= I 10 z dt (2.42) 



This model is considered to be an approximation of the Arrhenius model, which is 
valid over a wide range (4 to 160°C) of temperatures 133 : 



rl 



R 



( 



1 

T Z 



PE = —\e v /v 0J dt (2.43) 



where PE = integrated lethal effect, or pasteurization effect; E a = energy of activa- 
tion, J mol -1 ; R = 8.314 J mol -1 K _1 ; T= temperature, K; T = reference temperature, 
345 K; t = time, s; t = reference time, 15 s. The reference temperature (345 K or 
72°C) and time (15 s) correspond to the International Dairy Federation standard for 
pasteurization. 134 

It is often necessary for food processors to demonstrate that the process they 
wish to use is effective in delivering the required lethal effect for the product and 
microorganism of concern. The integrated lethal effect is a useful concept, as it 
allows two or more processes that use different time/temperature combinations to 
be compared for efficacy against food-borne pathogens; however, there are few data 
available relating microbial survival to processing conditions. This is of particular 
concern in the case of pasteurization of milk, where the only accepted test for proper 



2004 by Robin C. McKellar and Xuewen Lu 












1237_C02.fm Page 47 Wednesday, November 12, 2003 12:34 PM 







pasteurization is the alkaline phosphatase (AP) test. The relationship between AP 
inactivation and survival of food-borne pathogens is largely unknown, as is the 
response of AP to processing under alternative time/temperature combinations. 
Thus, modeling of HTST pasteurization of milk was studied extensively by McKel- 
lar and coworkers. 

The residence times in each section of a pilot-scale HTST pasteurizer and in 
each of six holding tubes with nominal holding times of 3 to 60 s were calibrated 
using the standard salt test. Temperatures were taken at the beginning and end of 
each section using thermocouples. The PE could then be determined for each selected 
holding time/temperature combination using Equation 2.43. Raw milk at a constant 
flow rate was allowed to equilibrate at each time/temperature, and a sample was 
taken at the outflow for analysis. Residual enzyme activity or microbial survivors 
were matched with the corresponding PE for fitting. 97 

The fitting was accomplished using an iterative procedure in which the log 10 % 
initial activity or viable counts were regressed on PE, with the value of EJR varied 
to minimize the error sum of squares. Nonlinearity (generally concavity) in the data 
was modeled using a power transformation (Equation 2.25). The final model was 
of the form 97 : 



log 10 % initial activity = a + b- PE ( 



(2.44) 




where a = intercept, b = slope, and c = power. Generally, the parameter estimates 
for at least three trials were pooled, and the model for AP is shown in Table 2.5. 

There is also a need to develop models for milk enzymes that might be used to 
confirm processing at temperatures above or below pasteurization. Higher temper- 
atures (>75°C) are appropriate for processing of more viscous products (such as ice- 
cream mix), while temperatures below pasteurization (63 to 65°C; termed subpas- 
teurization or thermization temperatures) are used to extend the storage life of bulk 




TABLE 2.5 

Model Parameters for Inactivation of Various Milk Enzymes 

and Food-Borne Pathogens during High-Temperature Short-Time 

Pasteurization 



Target 


Trials 


Intercept 


Slope 


Power 


EJR (x1 000) 


Alkaline phosphatase 


3 


2.05 


^.05 


0.50 


66.5 


y-Glutamyl transpeptidase 


3 


2.00 


-0.281 


0.75 


66.5 


Lactoperoxidase 


3 


2.12 


-0.10 


0.75 


59.0 


Catalase 


3 


1.94 


-2.65 


0.50 


82.0 


oc-L-Fucosidase 


3 


1.87 


-17.6 


1.00 


39.8 


Listeria innocua 


5 


1.86 


-24.9 


0.80 


59.5 


Listeria monocytogenes 


3 


1.68 


-18.4 


0.80 


48.5 


Enterobacter sakazakii 


3 


2.31 


-24.4 


0.65 


59.5 
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FIGURE 2.1 1 Linear model relating pasteurization effect (PE) and residual activity of lac- 
toperoxidase during high-temperature short-time (HTST) pasteurization of bovine milk. 

milk. Lactoperoxidase (LP) and y-glutamyl transpeptidase (TP) are two naturally 
occurring milk enzymes that are inactivated at higher temperatures. 135 Model param- 
eters for these two enzymes are given in Table 2.5. An example of an inactivation 
curve for LP is given in Figure 2.11, with the dotted lines representing the 95% 
confidence limits. There is close agreement among the three trials plotted, a char- 
acteristic common for all enzyme models. Models have also been developed for 
catalase 104 and oc-L-fucosidase (FC), 136 which are appropriate for subpasteurization 
temperatures (Table 2.5). 

Survival models for several food-borne pathogens have also been derived. List- 
eria innocua, a nonpathogen, is often used as a substitute for L. monocytogenes in 
situations (such as food processing environments) where it would be undesirable to 
introduce pathogens. 137 A model developed for L. innocua (Table 2.5) was shown 
to underpredict inactivation of L. monocytogenes', thus predictions are "fail-safe." 138 
Enterococcus faecium, a nonpathogen, is also used as a model organism for patho- 
gens, particularly in Europe. 139 The inactivation curve for this microorganism devi- 
ated strongly from linearity, and there were large intertrial variations. Thus, a random 
coefficient model using Equation 2.28 was used to fit the data. 96 The average In D- 
values for the two populations were 0.825 and 2.856. Models were also generated 
for Enterobacter sakazakii, an "emerging" pathogen found contaminating infant 
formula. 140 Model parameters compared with those for L. monocytogenes (Table 2.5) 
showed that E. sakazakii was more sensitive to pasteurization. 

Linear models for milk enzymes were characterized by limited intertrial vari- 
ability (Figure 2.1 1). This allowed validation of models using data from other trials 
that were not used in the construction of the models. In contrast, considerable 
variation was noted in experiments with microorganisms; thus a different approach 
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FIGURE 2.12 Probability densities for a-L-fucosidase, Listeria monocytogenes, and alkaline 
phosphatase generated from linear models for high- temperature short-time (HTST) pasteur- 
ization using Analytica®, with holding temperature and time of 66°C for 16 s. 

was taken. Model parameters were incorporated into risk analysis software (@RISK, 
Palisade Corporation, Newfield, NY) as normal distributions, with means taken from 
Table 2.5 and standard deviations taken from the intertrial variations. When simu- 
lations were performed (1500 iterations), outcomes (log reduction in this case) were 
expressed as distributions. 

Simulated log reductions were generated for AP, FC, and L. monocytogenes 
using a holding time of 65°C/15 s (thermization), and the probability density func- 
tions are shown in Figure 2.12. These conditions resulted in a narrow band of 
probabilities for AP, with greater predicted range for both FC and L. monocytogenes. 
AP is not completely inactivated, while FC (a potential indicator of thermization) 
experiences a >2 log reduction in most iterations. The mean log reduction of L. 
monocytogenes under these conditions is >3. 

Models that can predict the probability of achieving a desired level of safety are 
an important addition to risk assessment models, which are still largely qualitative 
and based primarily on expert opinion (see Chapter 6 for a more complete discussion 
on expert systems). The pasteurization models described above have been incorpo- 
rated into the risk analysis software Analytica® (Lumina Decision Systems, Los 
Gatos, CA), a commonly used software for building risk assessment models for the 
food industry. These models are now being incorporated into the USDA's Pathogen 
Modeling Program (available from http://www.arserrc.gov/mfs/pathogen.htm). 

2.2.6.2 Alternative Technologies 

Thermal treatment has been the traditional method for processing of many foods; 
however, with the increased consumer demand for fresh, less processed foods, new 
technologies have evolved. Some of these are based on temperature, such as micro- 
wave, radio frequency (RF), and ohmic heating, while others depend on other forms 
of microbial inactivation, such as high pressure (HP), pulse electric field (PEF), 
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pulsed or ultraviolet light, and ultrasound. In 1998, the U.S. Food and Drug Admin- 
istration commissioned the Institute of Food Technologists to provide scientific 
review and analysis of issues in food safety, food processing, and human health. 
The first of these reports, entitled "Kinetics of Microbial Inactivation for Alternative 
Food Processing Technologies," was released in 2000 141 and is available at the 
following web site: http://vm.cfsan.fda.gov/~comm/ift-toc.html. Since this report 
comprehensively reviews the scientific literature and makes recommendations for 
future research, it is beyond the scope of this chapter to reproduce this body of work. 
Instead, several key areas will be highlighted. 

Many novel thermal technologies base their antimicrobial effect on temperature; 
thus inactivation of microorganisms can be modeled using the traditional calculations 
for D- value and z-value (see earlier). Processes that depend on other mechanisms 
of inactivation such as HP and PEF require modified equations with different param- 
eters. For example, HP effects on microbial population can be modeled using a 
function similar to the traditional Z)-value 142 : 



log 



' D^ 



K D RJ 



(P-P») 



(2.45) 



•R 




where D R = the decimal reduction time at a reference pressure P R and z R is the 
pressure required for a 1-log reduction in D- value. An alternative model has been 
proposed by Weemaes et al. 143 : 




\n(k) = \n(k R )- 



r V{P-P^ 



v 



RT 



A 



(2.46) 



J 



where k R is the reaction rate constant, and P R the reference pressure, V the activation 
volume constant, P the pressure, and T A the absolute temperature. With PEF pro- 
cessing, a model describing the influence of the electric field intensity on reduction 
of microbial population can be described, which is similar to those used for thermal 
and pressure processing: 



log 



' D^ 



K D RJ 



(E-E R ) 



(2.47) 



'E 



where D R is the decimal reduction time at a reference field intensity E R , and the 
electric field coefficient z E is the increase in the electric field intensity E required to 
reduce the D-value by 1-log. An alternative model based on the Fermi equation was 
proposed by Peleg 144 : 






1 



E-E„ 



(2.48) 
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where E d is the electric field intensity when the microbial population has been 
reduced to 50%, and K is a coefficient based on the slope of the survivor curve. A 
similar model was proposed by Hulsheger et al. 145 : 



N 



f \ 
t 



K f cj 



E-E c 



(2.49) 




where t is the treatment time, t c the minimum treatment time for inactivation, E c the 
minimum field strength for inactivation, and K is a specific rate constant. This 
function is similar to Equation 2.48 except that it also accounts for exposure time 
at a given electric field intensity. 

The Institute of Food Technologies report has raised a number of relevant issues 
that would benefit from some discussion here. Kinetic parameters for microbial 
populations exposed to thermal treatment are well documented and provide a good 
basis upon which to develop models for alternative thermal processes. The nonther- 
mal models described above assume that microbial inactivation is a first-order 
reaction; however, as mentioned earlier, there is little direct evidence supporting this 
view. It will be necessary to further evaluate the adequacy of linear survival models, 
and to hopefully develop a universal model applicable to both thermal and non- 
thermal processes. In addition, experimental protocols have been found to be inad- 
equate to provide statistically reliable parameters for microbial reduction resulting 
from exposure to alternative technologies. This is particularly a problem with high 
pressure processing, where data are needed at different pressures with control of 
temperature and product. The inactivation mechanism for thermal destruction of 
microbes is generally well known, and evidence for additional independent mecha- 
nisms with processes such as ohmic heating is still lacking. Further work is needed 
to elucidate the mechanism of inhibition with alternative treatments such as PEF 
and HP, and to assess possible synergistic effects between alternative technologies 
and temperature. 

2.2.7 Injury/Repair Models 

Almost without exception, available models for microbial growth and death have 
been developed using fully viable, unstressed cells; thus the resulting models rep- 
resent the idealized scenario. It is well known that bacterial cells exposed to some 
form of sublethal stress require an adaptation or recovery period prior to growth; 
however, mathematical models do not incorporate the influence of stress. This was 
emphasized in a study designed to model the evolution of a log phase in L. mono- 
cytogenes, induced by acid, alkaline, and osmotic shocks. 146 When lag -phase cells 
(which are more sensitive to environmental stress than stationary-phase cells) were 
exposed to changes in pH or increased levels of NaCl, the subsequent generation 
times predicted by commercially available software were shorter than the observed 
experimental generation times. 

The physiological events that account for microbial injury and repair are poorly 
understood; thus there have been very few attempts to apply mathematical models 
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FIGURE 2.13 Change in populations of uninjured (O), injured (□), and dead (•) cells 
during sublethal heating according to Equation 2.50. 




to the phenomena of bacterial cell injury and resuscitation. The models that do exist 
are of two general types: those that aim to quantitate the extent of injury with 
increased exposure to stress, and those that attempt to predict the time required for 
repair and recovery of viability. 

Several attempts have been made to model the extension of the lag phase in 
response to stress. A model to describe the relationship between lag prior to growth 
and stress duration was proposed by Breand et al. 147 This model was developed to 
reflect the observation that the lag increased with increasing stress duration, and 
then decreased to a minimum lag at longer stress times. The empirical model 
described the influence of stress on the lag with a linear function, followed by a 
logistic decrease. Cheroutre-Vialette and Lebert 148 proposed the use of a recurrent 
neural network to model the changes in lag phase and growth rate experienced by 
L. monocytogenes exposed to osmotic and pH shock. Lambert and van der Ouderaa 149 
compared the relative ability of the Bioscreen (see Chapter 1) and viable counts to 
quantitate the inactivation of microorganisms by disinfection. They proposed a 
simple first-order inactivation reaction with accumulation of injured cells prior to 
complete loss of viability: 




fc 



fc, 



A { -> A 2 -^ P 



(2.50) 





where A x are the uninjured cells, A 2 the injured cells, and P are the dead cells. k x 
and k 2 are rate constants for injury and death, respectively. Populations of viable, 
injured, and dead cells were simulated based on the data of Lambert and van der 
Ouderaa, 149 and are shown in Figure 2.13. These responses were confirmed using 
image analysis of colony sizes on agar plates; viable and injured cells could be 
distinguished on the basis of size. Colony size was also used to quantitate cells of 
L. monocytogenes that had been injured by exposure to heat or starvation. 150 The 
colony size distribution was normal for uninjured cells, but demonstrated a right- 
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hand skew with injured cells. Percent sublethal injury could be related to colony 
area using a linear function. More recent studies using the flow cytometer to measure 
the distribution of the lag times of individual cells of Lactobacillus plantarum also 
showed a deviation from normality with heat-treated cells. The extreme value dis- 
tribution was found to be the best function for fitting both injured and uninjured cells: 



F(x) = 1- exp 



( 



v 



exp 



x — a 



b 



(2.51) 



-J 




where a and b are unknown parameters. 

There are few studies that aim to model the recovery of cells from injury. Injured 
cells can be differentiated on the basis of increased sensitivity to selective media 
(e.g., 5% NaCl), and it is thus possible to develop models to predict the time required 
for cells to repair damage due to stress. This process is complicated by our lack of 
information on the true nature of injury in bacterial cells, and the mechanism by 
which cells recover. The two -compartment kinetic model developed by Hills and 
Mackey to describe bacterial growth 4 was revised and extended to account for cell 
injury and resuscitation. 5 In the revised model, there are rate constants for injury 
(R) and resuscitation (R T ), and parameters to describe the decrease (a) and increase 
(b) of the injury and resuscitation curves. 5 This model was used to fit data from the 
resuscitation of L. monocytogenes after exposure to sublethal heat. 151 It was shown 
that resuscitation could best be described with a reduced model with the parameter 
for increasing rate of recovery (b T ) eliminated. A quadratic regression model was 
subsequently derived that expressed the lag as a function of temperature and the 
initial number of injured cells. 151 

2.2.8 Combined Growth/Death Models 

There have been a limited number of attempts to combine growth and death functions 
into single models. These are often simply combinations of functions such as the 
Gompertz or logistic with their mirror images. For example, a two-term model 
describing the behavior of Lactobacillus spp. during the ripening of fermented 
sausage incorporated a Gompertz function for both growth and death. 152 In a similar 
fashion, the logistic function and its mirror image, the Fermi equation, have been 
combined. 21 ' 101 The latter model has been expanded to include a proposed distribution 
of cell resistances to stress, resulting in a death model that varies in shape. 21 The 
Baranyi model for growth was also combined with its mirror image to describe 
growth and death for Brochothrix thermosphacta. 115 In this model, a smoothing 
function was included to account for the transition between growth and death phases. 
Other combined functions have used simple exponential growth and 
decline. 68 ' 153 ' 154 In one of these models, 68 the lag phase preceding growth was handled 
by a first-order step that represented spore germination, repair, or adaptation. Jones 
et al. 153 described the adaptation of cells to growth as a transition between cells in 
two states, immature and mature. This model reduced to a simple balance between 
growth and death, with variations in division and mortality rates being described by 
empirical functions. 153 
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It is questionable if expressing death as a mirror image of growth is valid. 
There is little direct evidence that the lag phases preceding growth or death are 
due to similar physiological phenomena, although a convincing theoretical argu- 
ment has been offered in support of this hypothesis. 65 It seems likely, however, 
that the stationary phase of growth and the "tailing" phase of inactivation are the 
result of different physiological processes. Models that address growth and death 
as different processes, and attempt to describe the response of microbes to their 
environment in terms of transitions between states, would seem to be the most 
useful for future development. 68 ' 153 
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